function [HHstates,wage_moments,transitions,phl,plh]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,wi_1,ff_2,wf_2,fi_2,wi_2,Wff_h,Wfi_h,Wfn_h,Wif_h,Wnf_h,Wii_h,Win_h,Wni_h,Wnn_h,Wff_l,Wfi_l,Wfn_l,Wif_l,Wnf_l,Wii_l,Win_l,Wni_l,Wnn_l)


df_1=Par1(1);
di_1=Par1(2);
lnf_1=Par1(3);
lni_1=Par1(4);
lfi_1=Par1(5);
lif_1=Par1(6);
df_2=Par1(7);
di_2=Par1(8);
lnf_2=Par1(9);
lni_2=Par1(10);
lfi_2=Par1(11);
lif_2=Par1(12);
p_1=Par1(13);
q_1=Par1(14);
p_2=Par1(15);
q_2=Par1(16);
nu_hl=Par1(28);
nu_lh=Par1(29);

alphaf=Par1(32);
betaf=Par1(33);
alphai=Par1(34);
betai=Par1(35);

% Data simulation (MSM)
N_households=15000; % # households
maxtime=300;   % # quarters
rng('default'); % set seed

% Pre-allocation
value=zeros(N_households,maxtime);
state=zeros(N_households,maxtime);
wage_draw1=ones(N_households,maxtime);
wage_draw2=ones(N_households,maxtime);

% Health shocks are exogenous and independent of LM conditions so can model health shocks structure upfront
Hstate = zeros(N_households,maxtime);
health_draws = rand(N_households,maxtime);
% Hstate (0=good health, 1=bad health)
for i = 1:N_households
    for j = 2:maxtime  % in time=1 households are born in good health (Hstate=0)
        if health_draws(i,j) < nu_hl & Hstate(i,j-1)==0 % prob of moving from good to bad health
            Hstate(i, j) = 1;
            elseif health_draws(i,j) >= nu_hl & Hstate(i,j-1)==0 % prob of staying in good health
            Hstate(i, j) = 0;
        elseif health_draws(i,j) < nu_lh & Hstate(i,j-1)==1 % prob of moving from bad to good health
            Hstate(i, j) = 0;
        elseif health_draws(i,j) >= nu_lh & Hstate(i,j-1)==1 % prob of staying in bad health
            Hstate(i, j) = 1;
        end
    end
end

for iter=1:N_households
    
    time=1; % households are born in the situation 9: "NN" and Healthy (h)
    scale = lni_1 + lni_2 + lnf_1 + lnf_2 + (1 - lni_1)*(1 - lni_2)*(1-lnf_1)*(1-lnf_2);
    randnumber=scale*rand;
    
    if randnumber < lnf_1
        %wage_draw1(iter,time)=gendist(ff_1',1,1);
        urn = random('Beta', alphaf, betaf);
        transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
        wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
        wage_draw2(iter,time)=1;
        value(iter,time)=Wfn_h(wage_aux)*(Wfn_h(wage_aux)>Wnn_h)+Wnn_h*(Wfn_h(wage_aux)<=Wnn_h);
        state(iter,time)=3*(Wfn_h(wage_aux)>Wnn_h)+9*(Wfn_h(wage_aux)<=Wnn_h);
        wage_draw1(iter,time) = wage_aux*(Wfn_h(wage_aux)>Wnn_h) + (Wfn_h(wage_aux)<=Wnn_h);
        
    elseif randnumber >=lnf_1 && randnumber<(lni_1+lnf_1)
        %wage_draw1(iter,time)=gendist(fi_1',1,1);
        urn = random ('Beta', alphai, betai);
        transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
        wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
        wage_draw2(iter,time)=1;
        value(iter,time)=Win_h(wage_aux)*(Win_h(wage_aux)>Wnn_h)+Wnn_h*(Win_h(wage_aux)<=Wnn_h);
        state(iter,time)=7*(Win_h(wage_aux)>Wnn_h)+9*(Win_h(wage_aux)<=Wnn_h);
        wage_draw1(iter,time) = wage_aux*(Win_h(wage_aux)>Wnn_h) + (Win_h(wage_aux)<=Wnn_h);
        
    elseif randnumber >=(lni_1+lnf_1) && randnumber<(lnf_2+lni_1+lnf_1)
        wage_draw1(iter,time)=1;
        %wage_draw2(iter,time)=gendist(ff_2',1,1);
        urn = random ('Beta', alphaf, betaf);
        transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
        wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
        value(iter,time)=Wnf_h(wage_aux)*(Wnf_h(wage_aux)>Wnn_h)+Wnn_h*(Wnf_h(wage_aux)<=Wnn_h);
        state(iter,time)=5*(Wnf_h(wage_aux)>Wnn_h)+9*(Wnf_h(wage_aux)<=Wnn_h);
        wage_draw2(iter,time) = wage_aux*(Wnf_h(wage_aux)>Wnn_h) + (Wnf_h(wage_aux)<=Wnn_h);
        
    elseif randnumber >=(lnf_2+lni_1+lnf_1) && randnumber<(lni_2+lnf_2+lni_1+lnf_1)
        wage_draw1(iter,time)=1;
        urn = random ('Beta', alphai, betai);
        transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
        wage_aux =find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
        %wage_draw2(iter,time)=gendist(fi_2',1,1);
        value(iter,time)=Wni_h(wage_aux)*(Wni_h(wage_aux)>=Wnn_h)+Wnn_h*(Wni_h(wage_aux)<Wnn_h);
        state(iter,time)=8*(Wni_h(wage_aux)>=Wnn_h)+9*(Wni_h(wage_aux)<Wnn_h);
        wage_draw2(iter,time) = wage_aux*(Wni_h(wage_aux)>=Wnn_h)+(Wni_h(wage_aux)<Wnn_h);
    else
        wage_draw1(iter,time)=1;
        wage_draw2(iter,time)=1;
        value(iter,time)=Wnn_h;
        state(iter,time)=9;
        
    end
    
    
    time=2;
    
    while time<=maxtime

        
        
        
        %% health states and value functions (1-9 for h and 10-18 for l)
        
        if Hstate(iter,time) == 0 % good health
            Wff = Wff_h; Wfi = Wfi_h; Wfn = Wfn_h; Wif = Wif_h; Wii = Wii_h; Win = Win_h; Wnf = Wnf_h; Wni = Wni_h; Wnn = Wnn_h;
        end
        
        
        if Hstate(iter,time) == 1 % bad health
            Wff = Wff_l; Wfi = Wfi_l; Wfn = Wfn_l; Wif = Wif_l; Wii = Wii_l; Win = Win_l; Wnf = Wnf_l; Wni = Wni_l; Wnn = Wnn_l;
        end
        
        
           %% FF in t-1
        if (state(iter,time-1) == 1 || state(iter,time-1) == 10)
            scale = df_1 + df_2 + lfi_1 + lfi_2 + (1 - df_1)*(1-df_2)*(1 - lfi_1)*(1 - lfi_2);
            randnumber=scale*rand;
            if  (randnumber < df_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wnf(wage_draw2(iter,time-1));
                state(iter,time)=5*(1-Hstate(iter,time))+13*Hstate(iter,time);
            end
            if ( randnumber >= df_1 && randnumber<(df_1+df_2))
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wfn(wage_draw1(iter,time-1));
                state(iter,time)=3*(1-Hstate(iter,time))+12*Hstate(iter,time);
            end
            if ( randnumber >= (df_1+df_2) && randnumber<(lfi_1+df_1+df_2))
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value_aux=Wif(wage_aux,wage_draw2(iter,time-1)); value_aux1 = Wff(wage_draw1(iter,time-1),wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=4*(value_aux>value_aux1)+ 1*(value_aux<=value_aux1) + 9*Hstate(iter,time);
                if (rem(state (iter,time) - state(iter,time-1),9) ==0)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            
            if (randnumber >= (lfi_1+df_1+df_2) && randnumber<(lfi_2+lfi_1+df_1+df_2))
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                %wage_draw2(iter,time)=gendist(fi_2',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
          
                value_aux=Wfi(wage_draw1(iter,time-1),wage_aux); value_aux1 = Wff(wage_draw1(iter,time-1),wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=2*(value_aux>value_aux1)+ 1*(value_aux<=value_aux1) + 9*Hstate(iter,time);             

                if (rem(state (iter,time) - state(iter,time-1),9)==0)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if (randnumber>=(lfi_2+lfi_1+df_1+df_2))
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time) = Wff(wage_draw1(iter,time-1),wage_draw2(iter,time-1));
                state(iter,time)=1*(1-Hstate(iter,time))+10*Hstate(iter,time);
            end
        end
        %% FI em t-1
        if (state(iter,time-1)==2 || state(iter,time-1) == 11)
            scale = df_1 + di_2 + lfi_1 + lif_2 + (1 - di_2)*(1 -df_1)*(1-lfi_1)*(1 - lif_2);
            randnumber=scale*rand;
            
            if (randnumber < df_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wni(wage_draw2(iter,time-1));
                state(iter,time)=8*(1-Hstate(iter,time))+17*Hstate(iter,time);
            end
            if (randnumber >= df_1 && randnumber<(df_1+di_2))
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wfn(wage_draw2(iter,time-1));
                state(iter,time)=3*(1-Hstate(iter,time))+12*Hstate(iter,time);
            end
            if randnumber >= (df_1+di_2) && randnumber<(lfi_1+df_1+di_2)
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
           
                value_aux=Wii(wage_aux,wage_draw2(iter,time-1)); value_aux1 = Wfi(wage_draw1(iter,time-1),wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=6*(value_aux>value_aux1)+ 2*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            
            if randnumber >=(lfi_1+df_1+di_2) && randnumber<(lif_2+lfi_1+df_1+di_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                %wage_draw2(iter,time)=gendist(ff_2',1,1);
                urn = random ('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                
                value_aux=Wff(wage_draw1(iter,time-1),wage_aux); value_aux1 = Wfi(wage_draw1(iter,time-1),wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=1*(value_aux>value_aux1)+ 2*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end

            end
            if randnumber>=(lif_2+lfi_1+df_1+di_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)= Wfi(wage_draw1(iter,time-1),wage_draw2(iter,time-1));
                state(iter,time)=2*(1-Hstate(iter,time))+11*Hstate(iter,time);
                
            end
        end
        
        %% FN em t-1
        if (state(iter,time-1)==3 || state(iter,time-1)==12)
            scale = df_1 + lfi_1 + lni_2 +lnf_2 + (1-df_1)*(1-lfi_1)*(1-lni_2)*(1-lnf_2);
            randnumber=scale*rand;
            
            if randnumber < df_1*(1-p_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wnn;
                state(iter,time)=9*(1-Hstate(iter,time))+18*Hstate(iter,time);
            end
            if randnumber >= df_1*(1-p_2) && randnumber < df_1
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                %wage_draw2(iter,time)=gendist(fi_2',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                value_aux=Wni(wage_aux); value_aux1 = Wnn;
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=8*(value_aux>value_aux1)+ 9*(value_aux<=value_aux1) + 9*Hstate(iter,time);
                if state (iter,time) == 9 || state(iter,time)==18
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >= df_1 && randnumber<(lfi_1+df_1)
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                
                value_aux=Win(wage_aux); value_aux1 = Wfn(wage_draw1(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=7*(value_aux>value_aux1)+ 3*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if state (iter,time) == 3 || state(iter,time)==12
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lfi_1+df_1) && randnumber<(lnf_2+lfi_1+df_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                % wage_draw2(iter,time)=gendist(ff_2',1,1);
                urn = random ('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                
                
                value_aux=Wff(wage_draw1(iter,time-1),wage_aux); value_aux1 = Wfn(wage_draw1(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=1*(value_aux>value_aux1)+ 3*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
                
  
            end
            if randnumber >= (lfi_1+df_1 + lnf_2) && randnumber<(lni_2+lnf_2+lfi_1+df_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                %wage_draw2(iter,time)=gendist(fi_2',1,1);
                
                value_aux=Wfi(wage_draw1(iter,time-1),wage_aux); value_aux1 = Wfn(wage_draw1(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=2*(value_aux>value_aux1)+ 3*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
                
        
            end
            if randnumber >= (lni_2+lnf_2+lfi_1+df_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wfn(wage_draw1(iter,time-1));
                state(iter,time)= 3*(1-Hstate(iter,time))+12*Hstate(iter,time);
            end
        end
        %% IF em t-1
        if (state(iter,time-1)==4 || state(iter,time-1)==13)
            scale = df_2 + di_1 + lfi_2 + lif_1 + (1 - di_1)*(1 -df_2)*(1-lfi_2)*(1 - lif_1);
            randnumber=scale*rand;
            
            if randnumber < di_1
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wnf(wage_draw2(iter,time-1));
                state(iter,time)=5*(1-Hstate(iter,time))+14*Hstate(iter,time);
            end
            if randnumber >= di_1 && randnumber<(df_2+di_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Win(wage_draw2(iter,time-1));
                state(iter,time)=7*(1-Hstate(iter,time))+16*Hstate(iter,time);
            end
            if randnumber >= (df_2+di_1) && randnumber<(lif_1+df_2+di_1)
                % wage_draw1(iter,time)=gendist(ff_1',1,1);
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                
                value_aux=Wff(wage_aux, wage_draw2(iter,time-1)); value_aux1 = Wif(wage_draw1(iter,time-1),wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=1*(value_aux>value_aux1)+ 4*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
           
             
            end
            if randnumber >=(lif_1+df_2+di_1) && randnumber<(lfi_2+lif_1+df_2+di_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                %wage_draw2(iter,time)=gendist(fi_2',1,1);
                
                value_aux=Wii(wage_draw1(iter,time-1),wage_aux); value_aux1 = Wif(wage_draw1(iter,time-1),wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=6*(value_aux>value_aux1)+ 4*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
               
            end
            if randnumber>=(lfi_2+lif_1+df_2+di_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wif(  wage_draw1(iter,time),   wage_draw2(iter,time));
                state(iter,time)=4*(1-Hstate(iter,time))+13*Hstate(iter,time);
            end
        end
        %% NF em t-1
        if (state(iter,time-1)==5 || state(iter,time-1)==14)
            scale = df_2 + lfi_2 + lnf_1 + lni_1 + (1 - lni_1)*(1 - lnf_1)*(1-lfi_2)*(1 - df_2);
            randnumber=scale*rand;
            
            if randnumber < df_2*(1-p_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wnn;
                state(iter,time)=9*(1-Hstate(iter,time))+18*Hstate(iter,time);
            end
            if randnumber >= df_2*(1-p_1) && randnumber < df_2
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                
                value_aux=Win(wage_aux); value_aux1 = Wnn;
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=7*(value_aux>value_aux1)+ 9*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                
                if state (iter,time) == 9 || state (iter,time) == 18
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= df_2 && randnumber<(lfi_2+df_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                %  wage_draw2(iter,time)=gendist(fi_2',1,1);
                
                value_aux=Wni(wage_aux); value_aux1 = Wnf(wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=8*(value_aux>value_aux1)+ 5*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
                
               
            end
            if randnumber >= (lfi_2+df_2) && randnumber<(lnf_1+lfi_2+df_2)
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                %wage_draw1(iter,time)=gendist(ff_1',1,1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                
                 value_aux=Wff(wage_aux, wage_draw2(iter,time-1)); value_aux1 = Wnf(wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=1*(value_aux>value_aux1)+ 5*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
                
              
            end
            if randnumber >= (lnf_1+lfi_2+df_2) && randnumber<(lni_1+lnf_1+lfi_2+df_2)
                % wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                
                value_aux=Wif(wage_aux, wage_draw2(iter,time-1)); value_aux1 = Wnf(wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=4*(value_aux>value_aux1)+ 5*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
                
            end
            if randnumber >= (lni_1+lnf_1+lfi_2+df_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wnf(wage_draw2(iter,time-1));
                state(iter,time)=5*(1-Hstate(iter,time))+14*Hstate(iter,time);
            end
        end
        %% II em t-1
        if (state(iter,time-1)==6 || state(iter,time-1)==15)
            scale = di_1 + di_2 + lif_1 + lif_2 + (1-di_1)*(1-di_2)*(1-lif_1)*(1-lif_2);
            randnumber=scale*rand;
            
            if randnumber < di_1
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wni(wage_draw2(iter,time-1));
                state(iter,time)=8*(1-Hstate(iter,time))+17*Hstate(iter,time);
            end
            if randnumber >= di_1 && randnumber<(di_1+di_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Win(wage_draw1(iter,time-1));
                state(iter,time)=7*(1-Hstate(iter,time))+16*Hstate(iter,time);
            end
            if randnumber >= (di_1+di_2) && randnumber<(lif_1+di_1+di_2)
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                % wage_draw1(iter,time)=gendist(ff_1',1,1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                
                value_aux=Wfi(wage_aux, wage_draw2(iter,time-1)); value_aux1 = Wii(wage_draw1(iter,time-1),wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=2*(value_aux>value_aux1)+ 6*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
                
               
            end
            if randnumber >= (lif_1+di_1+di_2) && randnumber<(lif_2+lif_1+di_1+di_2)
                
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                
                value_aux=Wif(wage_draw1(iter,time-1),wage_aux); value_aux1 = Wii(wage_draw1(iter,time-1),wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=4*(value_aux>value_aux1)+ 6*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                     wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            
            if randnumber>=(lif_2+lif_1+di_1+di_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wii(wage_draw1(iter,time-1),wage_draw2(iter,time-1));
                state(iter,time)=6*(1-Hstate(iter,time))+15*Hstate(iter,time);
            end
        end
        %% IN em t-1
        if (state(iter,time-1)==7 || state(iter,time-1)==16)
            scale = di_1 + lif_1 + lnf_2 + lni_2 + (1 - di_1)*(1 - lif_1)*(1 - lnf_2)*(1- lni_1);
            randnumber=scale*rand;
            
            if randnumber < di_1*(1-q_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wnn;
                state(iter,time)=9*(1-Hstate(iter,time))+18*Hstate(iter,time);
            end
            if randnumber >= di_1*(1-q_2) && randnumber < di_1
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                % wage_draw2(iter,time)=gendist(fi_2',1,1);
                
                value_aux=Wni(wage_aux); value_aux1 = Wnn;
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=8*(value_aux>value_aux1)+ 9*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                
                if state (iter,time) == 9 || state (iter,time) == 18
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >= di_1 && randnumber<(lif_1+di_1)
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                %wage_draw1(iter,time)=gendist(ff_1',1,1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                
                value_aux=Wfn(wage_aux); value_aux1 = Win(wage_draw1(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=3*(value_aux>value_aux1)+ 7*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lif_1+di_1) && randnumber<(lnf_2+lif_1+di_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random ('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                %wage_draw2(iter,time)=gendist(ff_2',1,1);
                
                value_aux=Wif(wage_draw1(iter,time-1),wage_aux); value_aux1 = Win(wage_draw1(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=4*(value_aux>value_aux1)+ 7*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lnf_2+lif_1+di_1) && randnumber<(lni_2+lnf_2+lif_1+di_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                %wage_draw2(iter,time)=gendist(fi_2',1,1);
                
                value_aux=Wii(wage_draw1(iter,time-1),wage_aux); value_aux1 = Win(wage_draw1(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=6*(value_aux>value_aux1)+ 7*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lni_2+lnf_2+lif_1+di_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Win(wage_draw1(iter,time-1));
                state(iter,time)=7*(1-Hstate(iter,time))+16*Hstate(iter,time);
            end
        end
        %% NI em t-1
        if (state(iter,time-1)==8 || state(iter,time-1)==17)
            scale = di_2 + lif_2 + lni_1 + lnf_1 + (1 - di_2)*(1 - lif_2)*(1 - lni_1)*(1 - lnf_1);
            randnumber=scale*rand;
            
            if randnumber < di_2*(1-q_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wnn;
                state(iter,time)=9*(1-Hstate(iter,time))+18*Hstate(iter,time);
            end
            if randnumber >= di_2*(1-q_1) && randnumber < di_2
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value_aux=Win(wage_aux); value_aux1 = Wnn;
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=7*(value_aux>value_aux1)+ 9*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                
                if state (iter,time) == 9 || state (iter,time) == 18
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= di_2 && randnumber<(lif_2+di_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                %wage_draw2(iter,time)=gendist(ff_2',1,1);
                urn = random ('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                value_aux=Wnf(wage_aux); value_aux1 = Wni(wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=5*(value_aux>value_aux1)+ 8*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
               
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lif_2+di_2) && randnumber<(lnf_1+lif_2+di_2)
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                % wage_draw1(iter,time)=gendist(ff_1',1,1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                
                value_aux=Wfi(wage_aux, wage_draw2(iter,time)); value_aux1 = Wni(wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=2*(value_aux>value_aux1)+ 8*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lnf_1+lif_2+di_2) && randnumber<(lni_1+lnf_1+lif_2+di_2)
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                
                value_aux=Wii(wage_aux, wage_draw2(iter,time)); value_aux1 = Wni(wage_draw2(iter,time-1));
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=6*(value_aux>value_aux1)+ 8*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lni_1+lnf_1+lif_2+di_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wni(wage_draw2(iter,time-1));
                state(iter,time)=8*(1-Hstate(iter,time))+17*Hstate(iter,time);
            end
        end
        
        %% NN in t-1
        if (state(iter,time-1) ==9 || state(iter,time-1)==18)
            scale = lnf_1 + lni_1 + lni_2 + lnf_2 + (1 - lni_2)*(1 - lni_1)*(1 - lnf_2)*(1 - lnf_1);
            randnumber=scale*rand;
            
            if randnumber < lnf_1
                %wage_draw1(iter,time)=gendist(ff_1',1,1);
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                wage_draw2(iter,time)=1;
                
                value_aux=Wfn(wage_aux); value_aux1 = Wnn;
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=3*(value_aux>value_aux1)+ 9*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >=lnf_1 && randnumber<(lni_1+lnf_1)
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=1;
                
                value_aux=Win(wage_aux); value_aux1 = Wnn;
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=7*(value_aux>value_aux1)+ 9*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >=(lni_1+lnf_1) && randnumber<(lnf_2+lni_1+lnf_1)
                wage_draw1(iter,time)=1;
                %wage_draw2(iter,time)=gendist(ff_2',1,1);
                urn = random ('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                
                value_aux=Wnf(wage_aux); value_aux1 = Wnn;
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=5*(value_aux>value_aux1)+ 9*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >=(lnf_2+lni_1+lnf_1) && randnumber<(lni_2+lnf_2+lni_1+lnf_1)
                wage_draw1(iter,time)=1;
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                %  wage_draw2(iter,time)=gendist(fi_2',1,1);
                
                value_aux=Wni(wage_aux); value_aux1 = Wnn;
                value(iter,time)=value_aux*(value_aux>value_aux1)+ value_aux1*(value_aux<=value_aux1);
                state(iter,time)=8*(value_aux>value_aux1)+ 9*(value_aux<=value_aux1) + 9*Hstate(iter,time);

                if (rem((state (iter,time) - state(iter,time-1)),9) == 0)
                
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber>=(lni_2+lnf_2+lni_1+lnf_1)
                value(iter,time)=Wnn;
                state(iter,time)=9*(1-Hstate(iter,time))+18*Hstate(iter,time);
            end
        end
        
        
        time=time+1;
        
    end
    
end

%% Moments calculation

% use only last 2 quarters of simulation for each household

N_householdsH = sum(state(:,end-1) < 10)
N_householdsL = sum(state(:,end-1) > 9)

FF_h=sum(state(:,end-1)==1)/N_householdsH;
FI_h=sum(state(:,end-1)==2)/N_householdsH;
FN_h=sum(state(:,end-1)==3)/N_householdsH;
IF_h=sum(state(:,end-1)==4)/N_householdsH;
NF_h=sum(state(:,end-1)==5)/N_householdsH;
II_h=sum(state(:,end-1)==6)/N_householdsH;
IN_h=sum(state(:,end-1)==7)/N_householdsH;
NI_h=sum(state(:,end-1)==8)/N_householdsH;
NN_h=sum(state(:,end-1)==9)/N_householdsH;

FF_l=sum(state(:,end-1)==10)/N_householdsL;
FI_l=sum(state(:,end-1)==11)/N_householdsL;
FN_l=sum(state(:,end-1)==12)/N_householdsL;
IF_l=sum(state(:,end-1)==13)/N_householdsL;
NF_l=sum(state(:,end-1)==14)/N_householdsL;
II_l=sum(state(:,end-1)==15)/N_householdsL;
IN_l=sum(state(:,end-1)==16)/N_householdsL;
NI_l=sum(state(:,end-1)==17)/N_householdsL;
NN_l=sum(state(:,end-1)==18)/N_householdsL;
FX_l=FF_l+FI_l+FN_l;
IX_l=IF_l+II_l+IN_l;

FF_agg=sum(state(:,end-1)==1|state(:,end-1)==10)/N_households;
FI_agg=sum(state(:,end-1)==2|state(:,end-1)==11)/N_households;
FN_agg=sum(state(:,end-1)==3|state(:,end-1)==12)/N_households;
IF_agg=sum(state(:,end-1)==4|state(:,end-1)==13)/N_households;
NF_agg=sum(state(:,end-1)==5|state(:,end-1)==14)/N_households;
II_agg=sum(state(:,end-1)==6|state(:,end-1)==15)/N_households;
IN_agg=sum(state(:,end-1)==7|state(:,end-1)==16)/N_households;
NI_agg=sum(state(:,end-1)==8|state(:,end-1)==17)/N_households;
NN_agg=sum(state(:,end-1)==9|state(:,end-1)==18)/N_households;


%HHstates_reduced=[FF_h;FI_h;FN_h;IF_h;NF_h;II_h;IN_h;NI_h;NN_h;FX_l;IX_l;NF_l;NI_l;NN_l] % use some aggreg moments for bad health
HHstates=[FF_h;FI_h;FN_h;IF_h;NF_h;II_h;IN_h;NI_h;NN_h;FF_agg;FI_agg;FN_agg;IF_agg;NF_agg;II_agg;IN_agg;NI_agg;NN_agg]

  
%% SALARIES MOMENTS
wage_formal_1=wf_1(wage_draw1(:,end-1));
wage_informal_1=wi_1(wage_draw1(:,end-1));
wage_formal_2=wf_2(wage_draw2(:,end-1));
wage_informal_2=wi_2(wage_draw2(:,end-1));

indicewf_1=state(:,end-1)<=3|state(:,end-1)==10|state(:,end-1)==11|state(:,end-1)==12;
wage_formal_1=sort(wage_formal_1(indicewf_1));
indicewi_1=find(state(:,end-1)==4|state(:,end-1)==6|state(:,end-1)==7|state(:,end-1)==13|state(:,end-1)==15|state(:,end-1)==16);
wage_informal_1=sort(wage_informal_1(indicewi_1));
indicewf_2=find(state(:,end-1)==1|state(:,end-1)==4|state(:,end-1)==5|state(:,end-1)==10|state(:,end-1)==13|state(:,end-1)==14);
wage_formal_2=sort(wage_formal_2(indicewf_2));
indicewi_2=find(state(:,end-1)==2|state(:,end-1)==6|state(:,end-1)==8| state(:,end-1)==11|state(:,end-1)==15|state(:,end-1)==17);
wage_informal_2=sort(wage_informal_2(indicewi_2));

sdwi_1 = std(wage_informal_1);
sdwi_2 = std(wage_informal_2);
sdwf_1 = std(wage_formal_1);
sdwf_2 = std(wage_formal_2);

meanwf_1=mean(wage_formal_1);

p10wf_1=prctile(wage_formal_1,10);
p25wf_1=prctile(wage_formal_1,25);
p50wf_1=prctile(wage_formal_1,50);
p75wf_1=prctile(wage_formal_1,75);
p90wf_1=prctile(wage_formal_1,90);



meanwf_2=mean(wage_formal_2);

p10wf_2=prctile(wage_formal_2,10);
p25wf_2=prctile(wage_formal_2,25);
p50wf_2=prctile(wage_formal_2,50);
p75wf_2=prctile(wage_formal_2,75);
p90wf_2=prctile(wage_formal_2,90);

meanwi_1=mean(wage_informal_1);

p10wi_1=prctile(wage_informal_1,10);
p25wi_1=prctile(wage_informal_1,25);
p50wi_1=prctile(wage_informal_1,50);
p75wi_1=prctile(wage_informal_1,75);
p90wi_1=prctile(wage_informal_1,90);


meanwi_2=mean(wage_informal_2);

p10wi_2=prctile(wage_informal_2,10);
p25wi_2=prctile(wage_informal_2,25);
p50wi_2=prctile(wage_informal_2,50);
p75wi_2=prctile(wage_informal_2,75);
p90wi_2=prctile(wage_informal_2,90);


%wage_moments=[meanwi_1; meanwf_1; meanwi_2 ; meanwf_2; sdwi_1; sdwf_1; sdwi_2; sdwf_2; p10wi_1 ; p10wf_1; p10wi_2 ; p10wf_2;...
    %p25wi_1; p25wf_1; p25wi_2 ; p25wf_2; p50wi_1; p50wf_1; p50wi_2 ; p50wf_2; p75wi_1; p75wf_1; p75wi_2 ; p75wf_2;...
    %p90wi_1; p90wf_1 ; p90wi_2 ; p90wf_2]; % wages
wage_moments=[meanwi_1; meanwf_1; meanwi_2; meanwf_2; sdwi_1; sdwf_1; sdwi_2; sdwf_2]; % only using mean and stdev


% conditional transitions
formal_1_time1=(state(:,end-1)==1|state(:,end-1)==2|state(:,end-1)==3|state(:,end-1)==10|state(:,end-1)==11|state(:,end-1)==12);
informal_1_time1=(state(:,end-1)==4|state(:,end-1)==6|state(:,end-1)==7|state(:,end-1)==13|state(:,end-1)==15|state(:,end-1)==16);
nonemp_1_time1=(state(:,end-1)==5|state(:,end-1)==8|state(:,end-1)==9|state(:,end-1)==14|state(:,end-1)==17|state(:,end-1)==18);
formal_2_time1=(state(:,end-1)==1|state(:,end-1)==4|state(:,end-1)==5|state(:,end-1)==10|state(:,end-1)==13|state(:,end-1)==14);
informal_2_time1=(state(:,end-1)==2|state(:,end-1)==6|state(:,end-1)==8|state(:,end-1)==11|state(:,end-1)==15|state(:,end-1)==17);
nonemp_2_time1=(state(:,end-1)==3|state(:,end-1)==7|state(:,end-1)==9|state(:,end-1)==12|state(:,end-1)==16|state(:,end-1)==18);

formal_1_time2=(state(:,end)==1|state(:,end)==2|state(:,end)==3|state(:,end)==10|state(:,end)==11|state(:,end)==12);
informal_1_time2=(state(:,end)==4|state(:,end)==6|state(:,end)==7|state(:,end)==13|state(:,end)==15|state(:,end)==16);
nonemp_1_time2=(state(:,end)==5|state(:,end)==8|state(:,end)==9|state(:,end)==14|state(:,end)==17|state(:,end)==18);
formal_2_time2=(state(:,end)==1|state(:,end)==4|state(:,end)==5|state(:,end)==10|state(:,end)==13|state(:,end)==14);
informal_2_time2=(state(:,end)==2|state(:,end)==6|state(:,end)==8|state(:,end)==11|state(:,end)==15|state(:,end)==17);
nonemp_2_time2=(state(:,end)==3|state(:,end)==7|state(:,end)==9|state(:,end)==12|state(:,end)==16|state(:,end)==18);

Dfn_1=sum(formal_1_time1==1 & nonemp_1_time2==1)/sum(formal_1_time1==1);
Din_1=sum(informal_1_time1==1 & nonemp_1_time2==1)/sum(informal_1_time1==1);
Dnf_1=sum(nonemp_1_time1==1 & formal_1_time2==1)/sum(nonemp_1_time1==1);
Dni_1=sum(nonemp_1_time1==1 & informal_1_time2==1)/sum(nonemp_1_time1==1);
Dfi_1=sum(formal_1_time1==1 & informal_1_time2==1)/sum(formal_1_time1==1);
Dif_1=sum(informal_1_time1==1 & formal_1_time2==1)/sum(informal_1_time1==1);
Dfn_2=sum(formal_2_time1==1 & nonemp_2_time2==1)/sum(formal_2_time1==1);
Din_2=sum(informal_2_time1==1 & nonemp_2_time2==1)/sum(informal_2_time1==1);
Dnf_2=sum(nonemp_2_time1==1 & formal_2_time2==1)/sum(nonemp_2_time1==1);
Dni_2=sum(nonemp_2_time1==1 & informal_2_time2==1)/sum(nonemp_2_time1==1);
Dfi_2=sum(formal_2_time1==1 & informal_2_time2==1)/sum(formal_2_time1==1);
Dif_2=sum(informal_2_time1==1 & formal_2_time2==1)/sum(informal_2_time1==1);
Dni_1_Dfn_2=sum(nonemp_1_time1==1 & informal_1_time2==1 & formal_2_time1==1 & nonemp_2_time2==1)/sum(nonemp_1_time1==1);
Dni_1_Din_2=sum(nonemp_1_time1==1 & informal_1_time2==1 & informal_2_time1==1 & nonemp_2_time2==1)/sum(nonemp_1_time1==1);
Dni_2_Dfn_1=sum(nonemp_2_time1==1 & informal_2_time2==1 & formal_1_time1==1 & nonemp_1_time2==1)/sum(nonemp_2_time1==1);
Dni_2_Din_1=sum(nonemp_2_time1==1 & informal_2_time2==1 & informal_1_time1==1 & nonemp_1_time2==1)/sum(nonemp_2_time1==1);

transitions=[Dfn_1;Din_1;Dnf_1;Dni_1;Dfi_1;Dif_1;Dfn_2;Din_2;Dnf_2;Dni_2;Dfi_2;Dif_2;Dni_1_Dfn_2;Dni_1_Din_2;Dni_2_Dfn_1;Dni_2_Din_1];

phl= sum(Hstate(:,end-1)==0 & Hstate(:,end)==1)/N_householdsH;

plh= sum(Hstate(:,end-1)==1 & Hstate(:,end)==0)/N_householdsL;


end





